Delay estimation from noisy time series 



Tom Ohira 

Sony Computer Science Laboratory 
3-14-13 Higashi-gotanda, Shinagawa, 
Tokyo 14-1, Japan 

Ryusuke Sawatari 
Department of Computer Science 
Keio University 
3-14-1 Hiyoshi, Yokohama 223, Japan 
(February 1, 2008) 

We propose here a method to estimate a delay from a time series taking advantage of analysis of 
random walks with delay. This method is applicable to a time series coming out of a system which 
is or can be approximated as a linear feedback system with delay and noise. We successfully test 
the method with a time series generated by discrete Langevin equation with delay. 



Estimation of delay from a noisy time series has attracted much attention. Especially, when the time series is 
chaotic, estimation of delay has a practical motivation: time-delayed coordinates are typically used to estimate 
fractal dimensions and Lyapunov exponents. There are series of works considering the subject from this viewpoint 
jl]-Q . Another viewpoint is to consider that a noisy time series consists of underlying deterministic dynamics with past 
influence and a noise term. Some statisticians have taken this stand and devised methods of analysis, for example, 
using the generalized Langevin equation |5|,|| and fluctuation dissipation theorem 0. In physiological fields, a more 
specialized case of the feedback delay associated with control system has attracted great deal of interest. A series of 
attempts has been made to estimate delay from physiological experimental data (see e.g., PHl^D- 

Against this background, we present here a method of estimating delay from a time series which is or is approximately 
generated by a noisy linear feedback system. We take advantage of analysis of random walks whose transition 
probability depends on the walker's position in a fixed interval past. Such random walks are termed delayed random 
walks and were proposed recently as a platform on which to study systems with both noise and delay . We will 
describe each step of the method in a transparent manner for implementation into computer algorithms. The method 
is tested to show its effectiveness on several test time series generated by discrete Langevin equation with delay jl2) . 

Let us first describe the delayed random walk, on whose analysis we base our method for delay estimation. We 
consider a random walk which takes a unit step in a unit time. The delayed random walk we start with is an extension 
of a position dependent random walk whose step toward the origin is more likely when no delay exists. Formally, it 
has the following definition: 

P{X t+1 = n;X t+1 - T = s) 

= g(n - 1, s - l)P(X t = n - 1; A t+1 _ r = s; X t _ r = s - 1) 
+ g(n - 1, s + l)P(X t = n - 1; X t+1 _ T = s; AV r = s + 1) 
+ f(n + 1, s - l)P(X t = n + 1; X t+ i- T = s; AV T = s - 1) 

+ f(n + l,s+ l)P(X f = n + 1; A t+1 _ T = s; AVr = s + 1), (1) 
f(x,y) = -(l + ax + (3y) 

9(x,y) = - (l-ax-py) (2) 

where the position of the walker at time t is X t , and P(X tl = ui;X t2 = 112) is the joint probability for the walker 
to be at Mi and 112 at time t\ and £2, respectively. f(x,y) and g(x,y) are transition probabilities to take a step to 
the negative and positive directions respectively. Hence the transition probability depends on both the current and 
the t steps past positions of the walker. We note that the above definition is approximate: we are assuming that 
the probability for the walker to be at positions which violate the condition < f(x, y) < 1 is negligible. This is an 
extension of the delayed random walk model discussed in |Dj ]. 

Following the similar argument as JHJ], we can derive a set of coupled equations which the stationary correlation 
function of the model obeys: 
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K(0) = (1 - 2a)K(0) + 1 - 2(3K(t) 

K(u) = (1 - a)K(u - I) - 0K(t -{u - I)), (1 < u < r) 

AT(zt) = (1 - a)AT(u - 1) - (3K((u - 1) - r), (u > r) (3) 

We can solve this set of equations iteratively for K(u) given a, (3 and r; examples are shown in Figure 1. We note 
that the oscillatory solutions are seen with sufficiently large r, but the shapes of the curves are different for the cases 
of a > p and /3 > a {L3|. 

The delayed random walk model presented here can be considered a model of a linear delayed feedback system 
with noise in probability space. This can be seen more transparently by considering its counterpart in physical space, 
which is given as follows in discrete time with white noise £t . 

X t+1 - X t = -aX t - pX t . T + (Zt 1 Z u )=8(t 1 -t 2 ). (4) 

If the system which generates a time series is or is approximated as a noisy linear feedback system, we can use the 
above set of equations to estimate the delay and other parameters. The basic idea is to "tune" each parameter so that 
the correlation function from the time series numerically matches the solution of equation (^). We can derive several 
concrete algorithms of different approaches based on this basic idea. In the following we present one such method 
which is simple with respect to both its concept and its implementation. The concrete steps are as follows. 

(1) As a prerequisite, we need to have a stationary noisy time series and some physical assumption that it is or is 
approximately generated by a linear feedback system with delay. (Some aspects of a time series such as whether it is 
chaotic or not can be checked by already known methods 0.) 

(2) Construct the auto-correlation function C(u) from the time series. If it is oscillating with some C(u) < 0, 
we can go to the next step. (If not oscillating and Civ) > 0, it is still possible to "tune" parameters in principle. 
However, as other methods may be more appropriate H, we do not consider this case here.) An example is shown in 
Figure 2. 

(3) From {C(u)}, we generate a "normalized set". Decide on the unit step size, and normalize the correlation 
function by the following requirement derived from (0). 



Hence, we generate 



A'(0) - A'(l) = i (5) 



We assume that with correctly estimated parameters, K(u) generated this way obeys (||). 

(4) Estimate delay r e around the "first zero" Tj of the correlation function; ti is the smallest number such that 
K{ ri ) «0. 

(5) With estimated r e , generate the following two sets of ordered pair L\ = {(yi(u), Zi(u))} and L 2 = 
{{V2{u),z 2 (u))} from Kiu) 

. . Kiu) . . K(u + 1) 

K[u) K{u) 

(6) Plot L\ and L 2 . We use the following relation derived from (||): 

z\{u) = (1 - a)yi(u) - (3, z 2 (v) = -py 2 (u) - (1 - a) (9) 

Thus, our assumption here is that if we have a correct estimate of r then both plots will be fitted with a linear 
function whose slope and intercept will give us a and /3. Example of these plots are shown in Figure 3. 

(7) Compute x 2 error for each plot, and define 

E(r e )=x\+xl (10) 

Our best estimate of r is the one which minimizes E(r e ) near n. Corresponding a and (3 is obtained as described in 
(6) (Figure 4). 
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We have tried this algorithm on several test time series data generated by equation (|4|) with various parameter 
ranges, and sample results are shown in Table 1. We found that choice of the number of correlation function points 
used, u max , occasionaly affect our results. We heuristically chose u max at a value up to which the graph of C{u) is 
rather clear, typically about two to four times Tj. The estimates are quite reasonable as shown here and typically 
better than "first zero" estimate Tj. 

We have described a method of estimating parameters from a time series produced by a noisy linear feedback 
system with a single stable point using analysis of delayed random walks. As mentioned, other algorithms based on 
the same basic idea of "tuning" in to the correlation function can be devised. A scheme of cross-checking estimated 
parameters from these different algorithms is currently being investigated [p3f. In several fields, models have been 
constructed which include the effects of time delays [p 
critical examination of these models with extension of including noise effects by comparisons with experimental time 
series, especially near the equilibrium state of the system. We are currently involved in the application of this and 
similar algorithms to experimental time series from biological systems, such as posture control data |l7j, which can 
physically be assumed to have a delayed feedback. 

The authors would like to thank Prof. M. Tokoro of Keio University and Sony CSL for providing an opportunity 
for this collaborative work. 



14-19]. The method presented here could possibly help in the 



[1] F. Takens, Lect. Notes Math. 898, 366 (1981). 

[2] P. Grassberger and I. Procaccia, Physica D 9, 189 (1983); A. Wolf, J. B. Swift, H. L. Swinney and J. Vastano, Physica D 
16, 285 (1985). 

[3] M. J. Biinner, M. Popp, Th. Mayer, A. Kittel, and J. Parisi Phys. Rev. E.54, R3082 (1996). 

[4] M. Sano and Y. Sawada, Phys. Rev. Lett. 55, 1082 (1985); J. -P. Eckmann and D. Ruelle, Rev. Mod. Phys. 57, 617 (1985). 

T. Tanaka, K. Aihara, and M. Taki, Phys. Rev. E.54, 2122 (1996). 
[5] R. Kubo, in 1965 Tokyo Summer Lectures in Theoretical Physics, Part I, Many-Body Theory (Benjamin, New York, 1966), 

pp.1-16; Rep. Prog. Phys.29, 255 (1966). 
[6] H. Mori, Prog. Theor. Phys. 33, 423 (1965). 
[7] See e.g., Y. Okabe, Amer. Math. Soc. Transl. 161, 19 (1994). 
[8] J. Milton and A. Longtin, Vision Res. 30, 515 (1990). 
[9] T. Ohira and J. G. Milton, Phys. Rev. E.52, 3277 (1995). 
[10] C. Eurich and J. G. Milton, (in press for Phys. Rev. E.) 

[11] T. Ohira, Sony Computer Science Laboratory Technical Report SCSL-TR-96-014, 1996 (in press for Phys. Rev. E). 
[12] U. Kiichler and B. Mensch, Stochastics and Stochastics Reports 40, 23 (1992). 
[13] Details will be shown elsewhere. S. Sawatari and T. Ohira, (in preparation). 

[14] M. C. Mackey and L. Glass, Science 197, 287 (1977); A. Longtin and J. Milton, Biol. Cybern. 61, 51 (1989); 
[15] C. M. Marcus and R. M. Westervelt, Phys. Rev. A 39, 347 (1989). 

[16] E. Villermaux, Phys. Rev. Lett. 75, 4618 (1995); K. Khrustova, G. Veser, and A. Mikhailov, Phys. Rev. Lett. 75, 3564 

(1995); K. Ikeda and M. Matsumoto, Physica D 29 223 (1987). 
[17] J. J. Collins and C. J. DeLuca, Exp. Brain Res. 95, 308 (1993); J. J. Collins and C. J. DeLuca, Phys. Rev. Lett. 73, 764 

(1994). 




3 



FIG. 1. Stationary correlation function K(u) as a function of steps u iteratively obtained as the solution of equation (3). 
The parameters are set as (a,f3,r) = (o)(0.1, 0.15, 10) and (6)(0.2, 0.1, 20). 
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FIG. 2. An example of time series (a) and associated correlation function C(u) generated from equation (4)(b). The 
parameters are set as (a,/3, t) = (0.2,0.1,20). 
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FIG. 3. Examples of plots of Li (a,b,c) and L2 (d,e,f), with various estimation of r e for the time series shown in Fig. 2 
with (a,(3,r) = (0.2,0.1,20). The estimates are r e = (o)(d)15, (6)(e)20, (c)(/)25. 
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FIG. 4. An example of plots of E(r e ) with various estimates around r; = 14 for the time series shown in Fig. 2 with 
(a, [3, t) = (0.2, 0.1, 20). (b) is with finer scale around the minimum of E. 
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TABLE I. Table of estimation results. For a time series with parameters of (a,fa.r), our method generates a i; fa and r e , 
where a;,/?; are estimated from graph of Li, i = 1, 2. 
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